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1 Introduction 

In this report, we present numerical experiments with two particular ABS algorithms: 

(1) The Huang or modified Huang algorithm, 

(2) The implicit QR algorithm, 

These methods are used for solving the following problem: 

(a) Finding the least-squares solution of an overdetermined linear system Ax = b, where 
A G R m > n , beR m ,xeR n and m > n. 

The considered algorithms belong to the scaled ABS class of algorithms for solving linear systems, 
that is defined by the following scheme: 

(A) Let x% G R n be arbitrary and H\ G R n ' n be nonsingular arbitrary. Set i = 1. 

(B) Compute the residual r, = Axi — b. If r, = 0, then stop, Xi solves the problem. Otherwise 
compute Si = HiA T Vi, where vi G R n is arbitrary, save that v±, . . . ,Vi are linearly inde- 
pendent. If S{ ^ 0, then go to (C). If Sj = and rfvi = 0, then set Xi+i = x,- n H i+ i = Hi 
and go to (F). If Sj = and rfvi 7^ 0, then stop, the system is incompatible. 



(C) Compute the search vector pi by 



Pi = Hfzi, 



where Zi G R n is arbitrary, save that zj HiA T Vi 7^ 0. 

(D) Update the estimate of the solution by 

%i-\-l VCi &iPii 

where the stepsize Q4 is given by 

oci = rfvi/pfA T Vi. 

(E) Update the Abaffian matrix by 

H l+1 =Hi- H l A T v i wf Hi/wf HiA T Vi , 
where Wi G R n is arbitrary, save that wj HiA T Vi 7^ 0. 
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(F) If i = m, then stop, Xi + \ solves the problem. Otherwise increment the index i by one and 
go to (B). 

From (E), it follows by induction that Hi + \A T Vi = and HK. 1 W{ = 0, where V{ = [v±, . . . , 
and Wi = [wi, . . . ,W{\. One can show that the implicit factorization V^APi = Li holds, where 
Pi = \pii ■ ■ ■ iPi] an d L{ is nonsingular lower triangular. Moreover the general solution of the 
scaled subsystem V? Ax = V^b can be expressed in the form 

x = Xi + H? q, (1) 

where q € R n is arbitrary (see § for the proof). 

The basic ABS class is the subclass of the scaled ABS class obtained by taking vi = e^, e» 
being the i-th unitary vector (i-th column of the unit matrix). In this case, residual ri = Ax{ — b 
need not be computed in (B), i-th element rfei = afxi — bi suffices. 

This report is organized as follows. In Section 2, a short description of individual algorithms 
is given. Section 3 contains some details concerning test matrices and numerical experiments 
and diuscusses the numerical results. The Appendix contains listing of all numerical results. 
For other numerical results on ABS methods see [1,11]. For a listing of the used ABS codes see 
[12]. 



2 The tested ABS algorithms for overdetermined linear systems 

In this report, we deal with two basic algorithms, one belonging to the basic ABS class and 
one belonging to the scaled ABS class. To simplify description, we will assume that A G R m,n , 
m < n, has full row rank so that ^ in (B). 

The Huang algorithm is obtained by the parameter choices Hi = I, Vi = ej, Z( = Oj, W{ = a^. 
Therefore 

Pi = Hidi (2) 

and 

H l+1 = Hi- Pipf/afpi, (3) 
From (Q) and (||), it follows by induction that pi € Range(Ai) and 

H i+l = I-P i D- 1 pT. (4) 

where Ai = [a\,...,ai], Pi = \pi,...,pi] and Di = diag{a^p\, . . . ,ajpi). Moreover Hi is 
symmetric, positive semidefinite and idempotent (it is the orthogonal projection matrix into 
Null(Ai_i)). Since the requirement pi £ Null(Ai-i) is crucial, we can improve orthogonality by 
iterative refinement Pi = Hjai, j > 1 (usually j = 2), obtaining the modified Huang algorithm. 

The Huang algorithm can be used for finding the minimum-norm solution to the compatible 
underdetermined system Ax = b, i.e. for minimizing s.t. Ax = b. To see this, we use the 
Lagrangian function and convexity of ||x||. Then x is a required solution if and only if x = A T u 
for some u € R m or x 6 Range{A T ). But 

in 

Xm+l = x l a iPi 
i=l 

by (D) and p. L € Range(Ai) C Range(A^), so that if x\ = 0, then x m +i € Range(A T ). 

A short description of two versions of the Huang and modified Huang algorithms follows. 
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Algorithm 1 

(Huang and modified Huang, formula ©). 
Set x\ = and ill = I. 
For i = 1 to m do 

Set j?j = -ffjOj (Huang) or 

Pi = Hi(Hidi) (modified Huang), 

di = afpi and x i+ i = Xi - ((afxi - bi)/di)pi. 

If i < m, then set J3i+i = Hi — pipj /di. 

end do 

Algorithm 2 

(Huang and modified Huang, formula @). 
Set x\ = and Po empty. 
For i = 1 to m do 

Set Pi = (I- Pi-iD^pTjai (Huang) or 

Pi = (I — Pi-iD^pT^I - Pi^Dr^PFjoi (modified Huang), 
di = afpi and x i+ i = Xi - ((afxi - bi)/di)pi. 
If i < m, then set Pi = [Pj_i,pj]. 
end do 

The implicit QR algorithm is obtained by the parameter choices Hi = I, Vi = Api, z% = e^, 
lOj = e^. Since Zi = ej and Wi = Ci, the matrices i^i and Pi have the same structure as that in 
the implicit LU algorithm. Using the implicit factorization property V^APi = Li, we can see 
that Vj^iVi = V^Api = so that the vectors Vj, j = 1, . . . ,i, are mutually orthogonal. 

The implicit QR algorithm can be used for finding the least-squares solution of an overde- 
termined linear system Ax = b, where A € R m ' n , b € R m , x E R n and m > n, which is obtained 
after at most n steps. To see this, we recall that algorithms from the scaled ABS class solve 
the scaled system Ax = V^b. Since V n = AP n , where P n is square nonsingular, the condi- 
tion Pj[A T Ax = P^A T b implies A T Ax = A T b, which defines the least-squares solution of the 
overdetermined system Ax = b. 

A short description of the implicit QR algorithm follows. 

Algorithm 3 

(Implicit QR). 

Set x\ = 0, r\ = —b and H\ = I, 
For i = 1 to m do 

Set pi = Hfa, Vi = Api 

(only (i — 1) (n — i + 1) nonzero elements of Hi is used), 
cq = rfvi/vfvi, x i+ i = Xi - ctiPi and r i+i = - a^. 
(only i nonzero elements of Xi are updated). 
If i < m, then set s» = Hf A T Vi and P^+i = Hi — SipJ /vjvi 
(only i{n — i) nonzero elements of Hi + \ are updated), 
end do 

We now consider the problem of solving overdetermined linear systems, in the least squares 
sense. Consider the linear system Ax = b, where A £ R m ' n , b £ R m , x € R n and m > n. This 
system can be solved in the least squares sense in n steps by the inplicit QR algorithm as was 
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shown above. Another possibility is based on the application of the Huang algorithm. Since the 
least-squares solution has to satisfy the normal equation 

A T Ax = A T b, 

it can be obtained from the solution of the following augmented system . 

Ax = y, (5) 

A T y = A T b. (6) 

Since y € Range(A) by (|5|), it is a minimum-norm solution of (||) and it can be obtained by the 
Huang algorithm. Having y, the compatible system (g) can be solved by any ABS algorithm. 
Moreover, using the implicit factorization property A T P n = L n , we can write 

L T n x = Pi Ax = Ply = b 

which as a system with a lower triangular matrix can be easily solved by the back substitution. 
The lower triangular matrix L n can be obtained as a by-product of the Huang algorithm. If we 
use formula (||), then di is i-th diagonal element of L n and P^a, contains the other z-th-column 
elements of L n (di is i-th column of the matrix A). 

Matrix L n has not to be stored since it can be reconstructed from columns of the matrix P n . 
However, the back substitution has to be realized in a slightly different way in this case. 

A short description of two versions of the Huang algorithm for least-squares solution of 
overdetermined systems follows. 

Algorithm 4 

(Huang and modified Huang for least-squares with stored L n ). 
Set x\ = and Po empty. 
For z = l to n do 

If i > 1, then compute gi = Pf_idi and copy it into the lower triangular matrix L n . 
Set pi = di - Pi-iD^gi (Huang) or 

Pi = (I - Pi-tDr^pT^ai - P^D^gi) (modified Huang). 

Set d{ = djpi and copy it into the lower triangular matrix L n . Set bi = b T pi. 

If i < n, then set Pi = [Pi-i,pi\. 

end do 

Solve the triangular system L^x = b. 
Algorithm 5 

(Huang and modified Huang for least-squares without stored L n ). 
Set x\ = and Po empty. 
For i = 1 to n do 

Set Pi = {I- Pi-iD^pT^di (Huang) or 

Pi = (I - Pi^D^pTjtf - Pi^Dr^Pl^di (modified Huang). 
Set di = djpi. 

If i < n, then set p = [Pj_i,pj]. 

end do 

Set /„ = b. 

For i = n to 1 do 

Set Xi = pffi/di. 

If i > 1, then set = fi — Xidi 

end do 
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3 The computational experiments 



Performance of ABS algorithms has been tested by using several types of ill-conditioned matrices. 
These matrices can be classified in the following way. The first letter distinguishes matrices with 
integer T and real 'R' elements, both actually stored as reals in double precision arithmetic. The 
second letter denotes randomly generated matrices 'R' or matrices determined by an explicit 
formula 'D'. For randomly generated matrices, a number specifying the interval for the random 
number generator follows, while the name of matrices determined by the explicit formula contains 
the formula number (Fl, F2, F3). The last letter of the name denotes a way for obtaining 
ill-conditioned matrices: 'R' - matrices with nearly dependent rows, 'C - matrices with nearly 
dependent columns, 'S' - nearly singular symmetric matrices, 'B' - both matrices in KKT system 
ill-conditioned. More specifically: 

IR500 Randomly generated matrices with integer elements uniformly distributed in the in- 
terval [-500,500]. 

IR500R Randomly generated matrices with integer elements uniformly distributed in the in- 
terval [-500,500] perturbed in addition to have two rows nearly dependent. 
IR500C Randomly generated matrices with integer elements uniformly distributed in the in- 
terval [-500,500] perturbed in addition to have two columns nearly dependent. 
RR100 Randomly generated matrices with real elements uniformly distributed in the interval. 
[-100,100] 



IDF1 
IDF2 
IDF3 
IR50 



Matrices with elements a 



Matrices with elements a 
Matrices with elements on 



\i — j\, 1 < i < m, 1 < j < n (Micchelli-Fiedler matrix). 
\i — j\ 2 , I < i < m, 1 < j < n. 
ij — \i + j — (m + n)/2\, 1, < i < m, 1 < j < n. 
Randomly generated matrices with integer elements uniformly distributed in the in- 
terval [-50,50]. 

Matrices with linearly dependent rows were obtained in the following way. The input data 
contain four integers which specify two row indices i±, 12, one column index 23 and one exponent 
14. Then the row is copied into aj 2 . Furthermore is set to zero and ai 2 i z to 2~ H . Similar 
procedures are used for columns and symmetric matrices. 

Solution vectors were generated randomly with integer or real elements uniformly distributed 
in the interval [-10,10]. Right hand sides for overdetermined systems were obtained by the 
following way. First, vector b was generated randomly with integer or real elements uniformly 
distributed in the interval [-10,10]. Then its first element together with elements in the first 
row of the matrix A were redefined by the formulas 61 = —1 and a\j = J27=2 a ijbj so that 
A T b = 0. The right-hand side vector was determined by the formula b = b + Ax*. Since 
A 1 " Ax* = A T (b — b) = A T b, the normal equation is satisfied and x* is a least squares solution of 
the system Ax = b. Vector b was generated nonzero, so that the system Ax = b is incompatible. 

We have tested the following particular algorithms: 

impl.qr5 The implicit QR algorithm: Algorithm 5 (subroutine alg5d.f). 

huang6 The Huang algorithm for overdetermined systems: Algorithm 6 (subroutine. 

alg6d.f). 

mod.huang6 The modified Huang algorithm for overdetermined systems: Algorithm 6 (sub- 
routine alg6d.f). 

huang7 The Huang algorithm for overdetermined systems: Algorithm 7 (subroutine 

alg7d.f). 

mod.huang7 The modified Huang algorithm for overdetermined systems: Algorithm 7 (sub- 
routine alg7d.f). 
qr lapack Subroutine DGELS from the LAPACK package. 
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svd lapack Subroutine DGELSS from the LAPACK package, 
gqr lapack Subroutine DGELSX from the LAPACK package. 

Notice that ABS algorithms were implemented in their basic form without partitioning into 
block or other special adjustements serving for speed increase as done in the LAPACK software. 

Detailed results of computational experiments are presented in the Appendix. For each 
selected problem, the type of matrix and the dimension is given. Furthermore, both the solu- 
tion and the residual errors together with the detected rank and the computational time are 
given for each tested algorithm. Computational experiments were performed on a Digital Unix 
Workstation in the double precision arithmetic (machine epsilon equal to about 10 -16 ). 

We have tested overdetermined systems with n » to/2, n = to/2 and n « to/2. These 
systems were solved by using Algorithms 5 - Algorithm 7 together with explicit QR and SVD 
decomposition based methods taken from the LAPACK package. 

The following tables give synthetic results for the 21 tested problems, the number at the 
intersection of the i-th row with the fc-th column indicating how many times the algorithm at 
the heading of the i-th row gave a lower error than the algorithm at the heading of the fc-th row 
(in case there is a second number, this counts the number of cases when difference was less that 
one percent). 



solution error - 21 overdetermined linear systems 



methods 


huang 


mod. 


huang 


mod. 


impl 


• qr 


svd 


gqr 




M 


6 


huang 


7 


huang 


qr5 


lap 


lap 


lap 




M 




6 




7 




ack 


ack 


ack 


total 


huang6 




4 


0/21 


4 


5 


9 


2 


2 


26/21 


mod . huang6 


17 




17 


1/12 


13/6 


16 


11 


8 


83/18 


huang7 


0/21 


4 




4 


5 


9 


2 


2 


26/21 


mod . huang7 


17 


8/12 


17 




14/6 


18 


11 


10 


95/18 


impl . qr5 


16 


2/6 


16 


1/6 




10 


8 


6 


59/12 


qr lapack 


12 


5 


12 


3 


11 




6/4 


9 


58/4 


svd lapack 


19 


10 


19 


10 


13 


11/4 




4/12 


86/16 


gqr lapack 


19 


13 


19 


11 


15 


12 


5/12 




94/12 
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residual error - 21 overdetermined linear systems 



methods 


huang 


mod. 


huang 


mod. 


impl . 


qr 


svd 


gq r 




ii 


6 


huang 


7 


huang 


qr5 


lap 


lap 


lap 








6 




7 




ack 


ack 


ack 


total 


huang6 




11 


4 


7 


17 


16 


9 


13 


77 


mod . huang6 


10 




6/2 


9 


15/1 


13 


10 


11/1 


74/4 


huang7 


17 


13/2 




9 


16 


16 


12 


12 


95/2 


mod . huang7 


14 


12 


12 




17 


14/1 


13 


13 


95/1 


impl . qr5 


4 


5/1 


5 


4 




8 


7 


7 


40/1 


qr lapack 


5 


8 


5 


6/1 


13 




9 


9 


55/1 


svd lapack 


12 


11 


9 


8 


14 


12 




13 


79 


gqr lapack 


8 


9/1 


9 


8 


14 


12 


8 




68/1 



Prom the results in the Appendix and the above tables we can state the following conclusions: 

(1) When solving well-conditioned overdetermined systems, the explicit QR algorithms based 
on the Householder orthogonalization are usually faster than ABS methods tested (es- 
pecially if n » to/2). Therefore, we cannot recommend the later for solving standard 
problems. 

(2) Interesting results were obtained for extremely ill-conditioned problems. The modified 
Huang and the implicit QR algorithms failed to solve systems with the matrix IDF2. 
Such systems was not solved by simple explicit QR methods as well. On the other hand, 
the modified Huang and the implicit QR algorithms found solutions of systems with the 
matrix IDF3 extremely fast and, moreover, they were able to determine the numerical 
rank correctly. Performance of these algorithms in that particular case is remarkable, they 
are at least 20 times faster then the best LAPACK routines. 

(3) The LAPACK routine DGELSS based on the SVD decomposition is extremely slow (es- 
pecially if n » to/2), not suitable for solving our problems (it can be substituted by the 
much faster routine DGELSX). 

(4) In term of solution error mod-huang7 and gqr are the best; in term of residual error huang7 
and mod.huang7 are the best. 
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4 Appendix: Test results on overdetermined linear systems 



matrix dimension method solution residual rank time 

m n error error 



IR500 


1050 


950 


huang6 


0, 


. 17D- 


■01 


0. 


,23D- 


14 


950 


32. 


.00 


IR500 


1050 


950 


mod . huang6 


0. 


78D- 


•10 


0. 


20D- 


15 


950 


63. 


.00 


IR500 


1050 


950 


huang7 


0. 


. 17D- 


•01 


0. 


16D- 


14 


950 


30. 


.00 


IR500 


1050 


950 


mod . huang7 


0, 


79D- 


•10 





.85D- 


15 


950 


60. 


.00 


IR500 


1050 


950 


impl . qr5 


0. 


51D- 


•09 


0, 


, 12D- 


13 


950 


57. 


.00 


IR500 


1050 


950 


qr lapack 


0. 


84D- 


09 





.83D- 


14 


950 


16 


.00 


IR500 


1050 


950 


svd lapack 


0. 


31D- 


■08 





.80D- 


14 


950 


119. 


.00 


IR500 


1050 


950 


gqr lapack 


0. 


,28D- 


■08 





.28D- 


14 


950 


26 


.00 


;ondition number: 


0.42D+09 




















IR500 


1400 


700 


huang6 


0, 


. 13D- 


•02 


0. 


.51D- 


14 


700 


23. 


.00 


IR500 


1400 


700 


mod . huang6 


0. 


41D- 


•11 


0. 


. 12D- 


13 


700 


49. 


.00 


IR500 


1400 


700 


huang7 


0. 


13D- 


•02 


0. 


,82D- 


14 


700 


24. 


.00 


IR500 


1400 


700 


mod . huang7 


0. 


39D- 


•11 


0. 


. 15D- 


14 


700 


47. 


.00 


IR500 


1400 


700 


impl . qr5 


0. 


33D- 


•10 


0. 


.31D- 


13 


700 


40. 


.00 


IR500 


1400 


700 


qr lapack 


0. 


54D- 


•10 





.56D- 


14 


700 


15. 


.00 


IR500 


1400 


700 


svd lapack 


0. 


,54D- 


■10 





. 18D- 


13 


700 


67. 


.00 


IR500 


1400 


700 


gqr lapack 


0. 


23D- 


■09 





.77D- 


14 


700 


23 


.00 



condition number: 0.65D+08 
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Test results for overdetcrmined linear systems - continued 



matrix dimension method solution residual rank time 

m n error error 



IR500 


2000 


400 


IR500 


2000 


400 


IR500 


2000 


400 


IR500 


2000 


400 


IR500 


2000 


400 


IR500 


2000 


400 


IR500 


2000 


400 


IR500 


2000 


400 


condition number: 


IR500C 


1050 


950 


IR500C 


1050 


950 


IR500C 


1050 


950 


IR500C 


1050 


950 


IR500C 


1050 


950 


IR500C 


1050 


950 


IR500C 


1050 


950 


IR500C 


1050 


950 


;ondition number: 


IR500C 


1400 


700 


IR500C 


1400 


700 


IR500C 


1400 


700 


IR500C 


1400 


700 


IR500C 


1400 


700 


IR500C 


1400 


700 


IR500C 


1400 


700 


IR500C 


1400 


700 


;ondition number: 


IR500C 


2000 


400 


IR500C 


2000 


400 


IR500C 


2000 


400 


IR500C 


2000 


400 


IR500C 


2000 


400 


IR500C 


2000 


400 


IR500C 


2000 


400 


IR500C 


2000 


400 



condition number: 



huang6 


0. 


35D-03 


mod . huang6 


0. 


96D-12 


huang7 


0. 


35D-03 


mod . huang7 


0, 


87D-12 


impl . qr5 


0, 


55D-11 


qr lapack 


0. 


. 18D-10 


svd lapack 


0. 


18D-10 


gqr lapack 


0. 


.65D-10 


. 24D+08 






huang6 


0. 


23D-01 


mod . huang6 


0. 


. 11D-01 


huang7 


0. 


23D-01 


mod . huang7 


0. 


11D-01 


impl . qr5 


0, 


12D-01 


qr lapack 


0. 


45D-02 


svd lapack 


0, 


69D+00 


gqr lapack 


0, 


69D+00 


0.25D+14 






huang6 


0. 


27D+00 


mod . huang6 


0. 


21D-05 


huang7 


0. 


27D+00 


mod . huang7 


0. 


21D-05 


impl . qr5 


0. 


, 12D+01 


qr lapack 


0. 


35D-06 


svd lapack 


0. 


35D-06 


gqr lapack 


0. 


. 11D-05 


0.51D+12 






huang6 


0. 


. 11D-01 


mod . huang6 


0. 


33D-03 


huang7 


0. 


. 11D-01 


mod . huang7 


0. 


33D-03 


impl . qr5 


0. 


80D-03 


qr lapack 


0, 


10D-02 


svd lapack 


0. 


.70D+00 


gqr lapack 


0, 


70D+00 



0.52D+13 



0, 


,95D- 


15 


400 


12. 


.00 





.78D- 


15 


400 


21. 


.00 


0. 


.41D- 


•15 


400 


11. 


.00 


0. 


.33D- 


15 


400 


21. 


.00 


0. 


. 11D- 


13 


400 


18. 


.00 


0. 


.46D- 


13 


400 


7. 


.00 


0. 


.54D- 


13 


400 


20. 


.00 


0, 


. 10D- 


12 


400 


12. 


.00 





. 15D- 


14 


950 


33. 


.00 


0. 


.19D- 


14 


950 


63. 


.00 


0. 


.36D- 


16 


950 


30. 


.00 


0. 


. 11D- 


•14 


950 


60. 


.00 


0. 


.40D- 
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Test results for overdetcrmined linear systems - continued 



matrix dimension method solution residual rank time 

m n error error 
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Test results for overdetcrmined linear systems - continued 



matrix dimension method solution residual rank time 

m n error error 
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Test results for overdetcrmined linear systems - continued 



matrix dimension method solution residual rank time 

m n error error 
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Test results for overdetcrmined linear systems - continued 



matrix dimension method solution residual rank time 

m n error error 
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